library(survival)
library(forestplot)
library(tidyverse)
library(survival)
library(regplot)
library(rms)
riskoas <- read.table("risk_riskscore_FPKM.txt",sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = 3)

fpkm <- read.table("clock_TCGALASSO.txt",sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = 1)


colnames(fpkm) <- substr(colnames(fpkm), 1, 16)
# 替换列名中的 "-" 为 "."
colnames(fpkm) <- gsub("-", ".", colnames(fpkm))
fpkm <- fpkm %>% t() %>% as.data.frame()
comgene <- intersect(rownames(fpkm),rownames(riskoas))
deres <- fpkm[comgene,]
a <- rownames(riskoas)
b <- rownames(deres)
identical(a,b)#对ginfo重排序
merged_data <- cbind(deres,riskoas)

##处理生存
library(survival)
library(forestplot)
library(tidyverse)
#下载生存信息
#xena官网：https://xenabrowser.net/datapages/?cohort=GDC%20TCGA%20Liver%20Cancer%20(LIHC)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443
#读取生存信息tsv文件
surv = read.table(file = 'TCGAclinical_all.txt', sep = '\t', header = TRUE) 
surv2 = read.table(file = 'TCGA-LGG.survival.tsv', sep = '\t', header = TRUE)
merged_df <- rbind(surv, surv2)
sorted_df <- merged_df[order(merged_df[,4], decreasing = FALSE),]
surv <- sorted_df

#整理生存信息数据
surv$sample <- gsub("-",".",surv$sample)
rownames(surv) <- surv$sample
surv <- surv[,-1]
surv <- surv[,-2]

comgene <- intersect(colnames(TCGACliLGG),colnames(TCGACliGBM))
TCGACliLGG1 <- TCGACliLGG[,comgene]
TCGACliGBM1 <- TCGACliGBM[,comgene]
comgene <- intersect(rownames(TCGACliLGG),rownames(surv))
a <- colnames(TCGACliLGG1)
b <- colnames(TCGACliGBM1)
identical(a,b)
merged_data2 <- rbind(TCGACliLGG1,TCGACliGBM1)

merged_data2$OS.time=merged_data2$OS.time/365   
rownames(merged_data2) <- substr(rownames(merged_data2), 1, 16)
clinical <- merged_data2[,c("days_to_last_follow_up", "vital_status", "gender", "race", "paper_IDH.status", "paper_Grade", "paper_X1p.19q.codeletion")]


write.table(cli,file="clinical_all.txt",sep="\t",quote=F,col.names=T)

cliFile="clinical_all.txt" 
cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)
cli=cli[apply(cli,1,function(x)any(is.na(match('unknow',x)))),,drop=F]
cli$age=as.numeric(cli$age)


risk <- merged_data

########
samSample=intersect(rownames(risk), rownames(cli))
risk1=risk[samSample,,drop=F]
cli=cli[samSample,,drop=F]
a <- rownames(risk1)
b <- rownames(cli)
identical(a,b)
rt=cbind(risk1, cli)
rt <- na.omit(rt)
write.table(rt,file="risk合并后TCGAclinical.txt",sep="\t",quote=F,col.names=T)
rt <- rt[,-6]
#????????ͼ
res.cox=coxph(Surv(OS.time, OS) ~ . , data = rt)
nom1=regplot(res.cox,
             plots = c("density", "boxes"),
             clickable=F,
             title="",
             points=TRUE,
             droplines=TRUE,
             observation=rt[16,],
             rank="sd",
             failtime = c(1,3,5),
             prfail = F)

#????ͼ???յ÷?
nomoRisk=predict(res.cox, data=rt, type="risk")
rt=cbind(risk1, Nomogram=nomoRisk)
outTab=rbind(ID=colnames(rt), rt)
write.table(outTab, file="nomoRisk.txt", sep="\t", col.names=F, quote=F)

#У׼????
pdf(file="calibration.pdf", width=5, height=5)
#1??У׼????
f <- cph(Surv(OS.time, OS) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=1)
cal <- calibrate(f, cmethod="KM", method="boot", u=1, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1),
     xlab="Nomogram-predicted OS (%)", ylab="Observed OS (%)", lwd=1.5, col="#33a02c", sub=F)
#3??У׼????
f <- cph(Surv(OS.time, OS) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=3)
cal <- calibrate(f, cmethod="KM", method="boot", u=3, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1), xlab="", ylab="", lwd=1.5, col="#3B4992FF", sub=F, add=T)
#5??У׼????
f <- cph(Surv(OS.time, OS) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=5)
cal <- calibrate(f, cmethod="KM", method="boot", u=5, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1), xlab="", ylab="",  lwd=1.5, col="#E64B35FF", sub=F, add=T)
legend('bottomright', c('1-year', '3-year', '5-year'),
       col=c("#33a02c","#3B4992FF","#E64B35FF"), lwd=1.5, bty = 'n')
dev.off()

